#include "IO2D.h"

#include <sstream>

#include "../Input/PicParams.h"
#include "../Parallel/ParallelMPI2D.h"
#include "../Field/Field2D.h"
#include "../PhysicalField/PhysicalField.h"
#include "../Species/Species.h"

using namespace std;

IO2D::IO2D(PicParams* params_, ParallelMPI* pmpi_) :
    IO(params_, pmpi_)
{

}

IO2D::~IO2D()
{
}


void IO2D::createFieldsGroup()
{
    fieldsGroup.dims_global[2] = params->dims_space_global[1] + 1;
    fieldsGroup.dims_global[1] = params->dims_space_global[0] + 1;
    fieldsGroup.dims_global[0] = 1;


    fieldsGroup.ndims_[0] = fieldsGroup.dims_global[0];
    fieldsGroup.ndims_[1] = fieldsGroup.dims_global[1];
    fieldsGroup.ndims_[2] = fieldsGroup.dims_global[2];

    fieldsGroup.offset[0] = 0;
    fieldsGroup.offset[1] = 0;
    fieldsGroup.offset[2] = 0;

    fieldsGroup.stride[0] = 1;
    fieldsGroup.stride[1] = 1;
    fieldsGroup.stride[2] = 1;

    fieldsGroup.block[0] = 1;
    fieldsGroup.block[1] = 1;
    fieldsGroup.block[2] = 1;

    //For attribute
    fieldsGroup.aDims = 3;


    //============ Create fieldsGroup ============================================

    addField(fields->rho_global);
    addField(fields->phi_global);
    addField(fields->Ex_global);
    addField(fields->Ey_global);
    addField(fields->rho_global_avg);
    addField(fields->phi_global_avg);
    addField(fields->Ex_global_avg);
    addField(fields->Ey_global_avg);
    for(int i = 0; i < fields->rho_s.size(); i++)
    {
        addField(fields->rho_s_global[i]);
        addField(fields->rho_s_global_avg[i]);
        addField(fields->Vx_s_global_avg[i]);
        addField(fields->Vy_s_global_avg[i]);
        addField(fields->Vz_s_global_avg[i]);
        addField(fields->Vp_s_global_avg[i]);
        addField(fields->T_s_global_avg[i]);

    }
    fieldsGroup.dataset_id.resize( fieldsGroup.dataset_stringName.size() );
}



void IO2D::createDiagsGroup()
{
    Diagnostic2D* diag2D = static_cast<Diagnostic2D*>(diag);

    string diag_name;
    const char* h5_name;
    hid_t dataset_id;
    int n_dim_data = 3;

    diagsGroup.dataset_stringName.resize(0);

    //=======set stride and block, and close dataset and group=================
    diagsGroup.stride[0] = 1;
    diagsGroup.stride[1] = 1;
    diagsGroup.stride[2] = 1;


    diagsGroup.block[0] = 1;
    diagsGroup.block[1] = 1;
    diagsGroup.block[2] = 1;

    //======= create diagsGroup ================================
    diag_name = "particleFlux";
    diagsGroup.dataset_stringName.push_back(diag_name);

    diag_name = "heatFlux";
    diagsGroup.dataset_stringName.push_back(diag_name);

    diag_name = "averageAngle";
    diagsGroup.dataset_stringName.push_back(diag_name);

    diag_name = "psiRate";
    diagsGroup.dataset_stringName.push_back(diag_name);
    
    diagsGroup.dataset_id.resize( diagsGroup.dataset_stringName.size() );
}



//write potential, rho and so on into hdf5 file every some timesteps
void IO2D::write( PicParams* params_, ParallelMPI* pmpi_, PhysicalField* fields, vector<Species*>& vecSpecies, Diagnostic* diag, int itime)
{
    const char* h5_name;
    int iDiag;
    int n_dim_data = 3;
    Diagnostic2D* diag2D = static_cast<Diagnostic2D*>(diag);


    if( itime % params->timestep_io == 0 && pmpi->isMaster() )
    {
        string step_output_string;

        step_output = itime / params->timestep_io;
        step_output_string = to_string(step_output);
        
        if(step_output_max >= 10 && step_output_max <100)
        {
            if(step_output < 10)
            {
                step_output_string = "0" + step_output_string;
            }
        }
        else if(step_output_max >= 100 && step_output_max <1000)
        {
            if(step_output < 10)
            {
                step_output_string = "00" + step_output_string;
            }
            else if(step_output<100)
            {
                step_output_string = "0" + step_output_string;
            }
        }
        else if(step_output_max >= 1000)
        {
            log_warning<<"step_output_max is too large, please change the code in IO2D.cpp";
        }

        data_file_name = "data/data" + step_output_string + ".h5";
        data_file_id = H5Fcreate( data_file_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);

        //============= write attributes, n_dim: number of dimension ======================
        hid_t attrs_dataspace_id, attrs_id;
        int n_dim = 2;
        hsize_t attrs_dims[1];
        attrs_dims[0] = 1;
        attrs_dataspace_id = H5Screate_simple(1, attrs_dims, NULL);
        attrs_id           = H5Acreate2(data_file_id, "n_dim", H5T_STD_I32BE, attrs_dataspace_id, H5P_DEFAULT, H5P_DEFAULT);
        H5Awrite(attrs_id, H5T_NATIVE_INT, &n_dim);
        H5Sclose(attrs_dataspace_id);
        H5Aclose(attrs_id);


        //=============write fields============================================
        fieldsGroup.group_id = H5Gcreate(data_file_id, "/Fields", H5P_DEFAULT, H5P_DEFAULT,H5P_DEFAULT);
        for(int i = 0; i < fieldsGroup.dataset_stringName.size(); i++)
        {
            fieldsGroup.dataspace_id = H5Screate_simple(n_dim_data, fieldsGroup.dims_global, NULL);
            h5_name = fieldsGroup.dataset_stringName[i].c_str();
            fieldsGroup.dataset_id[i] = H5Dcreate2(fieldsGroup.group_id, h5_name, H5T_NATIVE_DOUBLE, fieldsGroup.dataspace_id,
                                                    H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
            fieldsGroup.status = H5Dwrite(fieldsGroup.dataset_id[i], H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, fieldsGroup.dataset_data[i]);
            fieldsGroup.status = H5Sclose(fieldsGroup.dataspace_id);
            fieldsGroup.status = H5Dclose(fieldsGroup.dataset_id[i]);
        }
        fieldsGroup.status = H5Gclose(fieldsGroup.group_id);


        //=============write Diagnostics ============================================
        createDiagsGroup();

        diagsGroup.group_id = H5Gcreate(data_file_id, "/Diagnostic", H5P_DEFAULT, H5P_DEFAULT,H5P_DEFAULT);
        //particle flux
        iDiag = 0;
        diagsGroup.dims_global[0] = 1;
        diagsGroup.dims_global[1] = diag2D->n_species;
        diagsGroup.dims_global[2] = diag2D->n_segment2d_total;
        h5_name = diagsGroup.dataset_stringName[iDiag].c_str();
        diagsGroup.dataspace_id = H5Screate_simple(n_dim_data, diagsGroup.dims_global, NULL);
        diagsGroup.dataset_id[iDiag] = H5Dcreate2(diagsGroup.group_id, h5_name, H5T_NATIVE_DOUBLE, diagsGroup.dataspace_id,
                                                    H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        for(int iSpec = 0; iSpec < diag2D->n_species; iSpec++)
        {
            int offset_line = 0;
            for(int iLine = 0; iLine < diag2D->n_line; iLine++)
            {
                diagsGroup.offset[0] = 0;
                diagsGroup.offset[1] = iSpec;
                diagsGroup.offset[2] = offset_line;
                diagsGroup.count[0] = 1;
                diagsGroup.count[1] = 1;
                diagsGroup.count[2] = diag2D->n_segment2d[iLine];
                diagsGroup.memspace_id = H5Screate_simple(n_dim_data, diagsGroup.count, NULL);
                diagsGroup.status = H5Sselect_hyperslab(diagsGroup.dataspace_id, H5S_SELECT_SET, diagsGroup.offset,
                                                    diagsGroup.stride, diagsGroup.count, diagsGroup.block);
                diagsGroup.status = H5Dwrite(diagsGroup.dataset_id[iDiag], H5T_NATIVE_DOUBLE, diagsGroup.memspace_id,
                                        diagsGroup.dataspace_id, H5P_DEFAULT, &(diag2D->particleFlux_global[iSpec][iLine][0]));
                diagsGroup.status = H5Sclose(diagsGroup.memspace_id);
                offset_line += diag2D->n_segment2d[iLine];
            }
        }
        diagsGroup.status = H5Sclose(diagsGroup.dataspace_id);
        diagsGroup.status = H5Dclose(diagsGroup.dataset_id[iDiag]);

        //heat flux
        iDiag++;
        diagsGroup.dataspace_id = H5Screate_simple(n_dim_data, diagsGroup.dims_global, NULL);
        h5_name = diagsGroup.dataset_stringName[iDiag].c_str();
        diagsGroup.dataset_id[iDiag] = H5Dcreate2(diagsGroup.group_id, h5_name, H5T_NATIVE_DOUBLE, diagsGroup.dataspace_id,
                                                    H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        for(int iSpec = 0; iSpec < diag2D->n_species; iSpec++)
        {
            int offset_line = 0;
            for(int iLine = 0; iLine < diag2D->n_line; iLine++)
            {
                diagsGroup.offset[0] = 0;
                diagsGroup.offset[1] = iSpec;
                diagsGroup.offset[2] = offset_line;
                diagsGroup.count[0] = 1;
                diagsGroup.count[1] = 1;
                diagsGroup.count[2] = diag2D->n_segment2d[iLine];
                diagsGroup.memspace_id = H5Screate_simple(n_dim_data, diagsGroup.count, NULL);
                diagsGroup.status = H5Sselect_hyperslab(diagsGroup.dataspace_id, H5S_SELECT_SET, diagsGroup.offset,
                                                    diagsGroup.stride, diagsGroup.count, diagsGroup.block);
                diagsGroup.status = H5Dwrite(diagsGroup.dataset_id[iDiag], H5T_NATIVE_DOUBLE, diagsGroup.memspace_id,
                                        diagsGroup.dataspace_id, H5P_DEFAULT, &(diag2D->heatFlux_global[iSpec][iLine][0]));
                diagsGroup.status = H5Sclose(diagsGroup.memspace_id);
                offset_line += diag2D->n_segment2d[iLine];
            }
        }
        diagsGroup.status = H5Sclose(diagsGroup.dataspace_id);
        diagsGroup.status = H5Dclose(diagsGroup.dataset_id[iDiag]);

        //average angle
        iDiag++;
        diagsGroup.dataspace_id = H5Screate_simple(n_dim_data, diagsGroup.dims_global, NULL);
        h5_name = diagsGroup.dataset_stringName[iDiag].c_str();
        diagsGroup.dataset_id[iDiag] = H5Dcreate2(diagsGroup.group_id, h5_name, H5T_NATIVE_DOUBLE, diagsGroup.dataspace_id,
                                                    H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        for(int iSpec = 0; iSpec < diag2D->n_species; iSpec++)
        {
            int offset_line = 0;
            for(int iLine = 0; iLine < diag2D->n_line; iLine++)
            {
                diagsGroup.offset[0] = 0;
                diagsGroup.offset[1] = iSpec;
                diagsGroup.offset[2] = offset_line;
                diagsGroup.count[0] = 1;
                diagsGroup.count[1] = 1;
                diagsGroup.count[2] = diag2D->n_segment2d[iLine];
                diagsGroup.memspace_id = H5Screate_simple(n_dim_data, diagsGroup.count, NULL);
                diagsGroup.status = H5Sselect_hyperslab(diagsGroup.dataspace_id, H5S_SELECT_SET, diagsGroup.offset,
                                                    diagsGroup.stride, diagsGroup.count, diagsGroup.block);
                diagsGroup.status = H5Dwrite(diagsGroup.dataset_id[iDiag], H5T_NATIVE_DOUBLE, diagsGroup.memspace_id,
                                        diagsGroup.dataspace_id, H5P_DEFAULT, &(diag2D->averageAngle_global[iSpec][iLine][0]));
                diagsGroup.status = H5Sclose(diagsGroup.memspace_id);
                offset_line += diag2D->n_segment2d[iLine];
            }
        }
        diagsGroup.status = H5Sclose(diagsGroup.dataspace_id);
        diagsGroup.status = H5Dclose(diagsGroup.dataset_id[iDiag]);

        /*
        //psiRate
        iDiag++;
        diagsGroup.dataspace_id = H5Screate_simple(n_dim_data, diagsGroup.dims_global, NULL);
        h5_name = diagsGroup.dataset_stringName[iDiag].c_str();
        diagsGroup.dataset_id[iDiag] = H5Dcreate2(diagsGroup.group_id, h5_name, H5T_NATIVE_DOUBLE, diagsGroup.dataspace_id,
                                                    H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        for(int iPsi = 0; iPsi < diag2D->psiRate_global.size(); iPsi++)
        {
            diagsGroup.offset[0] = iPsi;
            diagsGroup.offset[1] = 0;
            diagsGroup.offset[2] = 0;
            diagsGroup.memspace_id = H5Screate_simple(n_dim_data, diagsGroup.count, NULL);
            diagsGroup.status = H5Sselect_hyperslab(diagsGroup.dataspace_id, H5S_SELECT_SET, diagsGroup.offset,
                                                diagsGroup.stride, diagsGroup.count, diagsGroup.block);
            diagsGroup.status = H5Dwrite(diagsGroup.dataset_id[iDiag], H5T_NATIVE_DOUBLE, diagsGroup.memspace_id,
                                    diagsGroup.dataspace_id, H5P_DEFAULT, (diag2D->psiRate_global[iPsi])->data_);
            diagsGroup.status = H5Sclose(diagsGroup.memspace_id);
        }
        diagsGroup.status = H5Sclose(diagsGroup.dataspace_id);
        diagsGroup.status = H5Dclose(diagsGroup.dataset_id[iDiag]);
        */

        diagsGroup.status = H5Gclose(diagsGroup.group_id);

        /*
        //write particle velocity distribution function
        ptclsGroup.group_id = H5Gcreate(data_file_id, "/VDF", H5P_DEFAULT, H5P_DEFAULT,H5P_DEFAULT);
        for(int i = 0; i < ptclsGroup.dataset_id.size(); i++)
        {
            ptclsGroup.dataspace_id = H5Screate_simple(n_dim_data, ptclsGroup.dims_global, NULL);
            h5_name = ptclsGroup.dataset_stringName[i].c_str();
            ptclsGroup.dataset_id[i] = H5Dcreate2(ptclsGroup.group_id, h5_name, H5T_NATIVE_DOUBLE, ptclsGroup.dataspace_id,
                                                    H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
            ptclsGroup.status = H5Dwrite(ptclsGroup.dataset_id[i], H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, ptclsGroup.dataset_data[i]);
            ptclsGroup.status = H5Sclose(ptclsGroup.dataspace_id);
            ptclsGroup.status = H5Dclose(ptclsGroup.dataset_id[i]);
        }


        ptclsGroup.status = H5Gclose( ptclsGroup.group_id );
        */


        status = H5Fclose(data_file_id);
    }



}

//read electric fields from hdf5 file for initial condition, or constant external field
void IO2D::read( PicParams* params_, ParallelMPI* pmpi_, PhysicalField* fields, vector<Species*>& vecSpecies, Diagnostic* diag, int itime)
{
    hid_t group_id, dataset_id;
    herr_t status;
    const char* h5_name;
    int iDiag;
    int n_dim_data = 3;

    if(pmpi->isMaster() )
    {
        data_file_name = params->read_data_file_name;
        data_file_id = H5Fopen( data_file_name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);

        group_id = H5Gopen(data_file_id, "/Fields", H5P_DEFAULT);

        //=============read fields phi ============================================
        dataset_id = H5Dopen(group_id, "phi_global_avg", H5P_DEFAULT);
        status = H5Dread (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, fields->phi_global);

        //=============read fields Ex ============================================
        dataset_id = H5Dopen(group_id, "Ex_global_avg", H5P_DEFAULT);
        status = H5Dread (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, fields->Ex_global);
        status = H5Dclose( dataset_id );

        //=============read fields Ey ============================================
        dataset_id = H5Dopen(group_id, "Ey_global_avg", H5P_DEFAULT);
        status = H5Dread (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, fields->Ey_global);
        status = H5Dclose( dataset_id );

        status = H5Gclose(group_id);
        status = H5Fclose(data_file_id);
    }
    pmpi->barrier();

    ParallelMPI2D* pmpi2D = static_cast<ParallelMPI2D*>(pmpi);

    Field2D<double>* Ex2D = static_cast<Field2D<double>*>(fields->Ex_);
    Field2D<double>* Ey2D = static_cast<Field2D<double>*>(fields->Ey_);
    Field2D<double>* phi2D = static_cast<Field2D<double>*>(fields->phi_);

    Field2D<double>* Ex2D_global    = static_cast<Field2D<double>*>(fields->Ex_global);
    Field2D<double>* Ey2D_global    = static_cast<Field2D<double>*>(fields->Ey_global);
    Field2D<double>* phi2D_global    = static_cast<Field2D<double>*>(fields->phi_global);

    pmpi2D->scatterField(Ex2D_global, Ex2D);
    pmpi2D->scatterField(Ex2D_global, Ex2D);
    pmpi2D->scatterField(phi2D_global, phi2D);
}

void IO2D::writeGrid(Grid* grid)
{
    hid_t       grid_dataspace_id;
    hid_t       grid_dataset_id;
    herr_t      grid_status;
    string      grid_dataset_name;

    hsize_t     count[3];              /* size of subset in the file */
    hsize_t     offset[3];             /* subset offset in the file */
    hsize_t     stride[3];
    hsize_t     block[3];

    int ii;
    int grid_ndim = 3;
    hsize_t grid_dims_global[3];

    Grid2D* grid2D = static_cast<Grid2D*>(grid);
    grid_file_name  = "data/grid.h5";
    grid_file_id    = H5Fcreate( grid_file_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);

    grid_dims_global[0] = 1;
    grid_dims_global[1] = grid2D->globalDims_[0];
    grid_dims_global[2] = grid2D->globalDims_[1];


    //=============write grid============================================
    grid_dataspace_id = H5Screate_simple(grid_ndim, grid_dims_global, NULL);
    grid_dataset_name = "is_wall";
    grid_dataset_id = H5Dcreate2(grid_file_id, grid_dataset_name.c_str(), H5T_NATIVE_INT, grid_dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    grid_status = H5Dwrite(grid_dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, grid2D->iswall_global_);
    grid_status = H5Sclose(grid_dataspace_id);
    grid_status = H5Dclose(grid_dataset_id);

    grid_dataspace_id = H5Screate_simple(grid_ndim, grid_dims_global, NULL);
    grid_dataset_name = "bndr_type";
    grid_dataset_id = H5Dcreate2(grid_file_id, grid_dataset_name.c_str(), H5T_NATIVE_INT, grid_dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    grid_status = H5Dwrite(grid_dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, grid2D->bndr_global_);
    grid_status = H5Sclose(grid_dataspace_id);
    grid_status = H5Dclose(grid_dataset_id);

    grid_dataspace_id = H5Screate_simple(grid_ndim, grid_dims_global, NULL);
    grid_dataset_name = "bndr_val";
    grid_dataset_id = H5Dcreate2(grid_file_id, grid_dataset_name.c_str(), H5T_NATIVE_DOUBLE, grid_dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    grid_status = H5Dwrite(grid_dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, grid2D->bndrVal_global_);
    grid_status = H5Sclose(grid_dataspace_id);
    grid_status = H5Dclose(grid_dataset_id);

    /*
    for(int iLine = 0; iLine < grid2D->lines.size(); iLine++)
    {
        for(int isegment2d = 0; isegment2d < grid2D->lines[iLine].size(); isegment2d++)
        {
            cout<<grid2D->lines[iLine][isegment2d].grid_point[0]<<" "<<grid2D->lines[iLine][isegment2d].grid_point[1]<<endl;
            cout<<grid2D->lines[iLine][isegment2d].start_point[0]<<" "<<grid2D->lines[iLine][isegment2d].start_point[1]<<endl;
            cout<<grid2D->lines[iLine][isegment2d].end_point[0]<<" "<<grid2D->lines[iLine][isegment2d].end_point[1]<<endl;
        }
        cout<<endl;
        cout<<"================================"<<endl;
    }
    */

    //===================== write boundary lines ======================================
    hid_t segment2d_type, memtype, point_type, grid_point_type, normal_type;
    hid_t segment2d_dataspace_id, segment2d_dataset_id;
    int ndim_segment2d;
    hsize_t dim1d_segment2d[1];
    hsize_t dim2d_segment2d[2];

    int *n_segment2d = new int[grid2D->n_line];
    for(int i = 0; i < grid2D->n_line; i++)
    {
        n_segment2d[i] = grid2D->n_segment2d[i];
    }

    double *start_point = new double[grid2D->n_segment2d_total * 2];
    ii = 0;
    for(int iLine = 0; iLine < grid2D->lines.size(); iLine++)
    {
        for(int isegment2d = 0;  isegment2d < grid2D->lines[iLine].size(); isegment2d++)
        {
            start_point[ii * 2]     = grid2D->lines[iLine][isegment2d].start_point[0];
            start_point[ii * 2 + 1] = grid2D->lines[iLine][isegment2d].start_point[1];
            ii++;
        }
    }

    double *end_point = new double[grid2D->n_segment2d_total * 2];
    ii = 0;
    for(int iLine = 0; iLine < grid2D->lines.size(); iLine++)
    {
        for(int isegment2d = 0;  isegment2d < grid2D->lines[iLine].size(); isegment2d++)
        {
            end_point[ii * 2]     = grid2D->lines[iLine][isegment2d].end_point[0];
            end_point[ii * 2 + 1] = grid2D->lines[iLine][isegment2d].end_point[1];
            ii++;
        }
    }

    double *length = new double[grid2D->n_segment2d_total];
    ii = 0;
    for(int iLine = 0; iLine < grid2D->lines.size(); iLine++)
    {
        for(int isegment2d = 0;  isegment2d < grid2D->lines[iLine].size(); isegment2d++)
        {
            length[ii] = grid2D->lines[iLine][isegment2d].length;
            ii++;
        }
    }

    int *grid_point0 = new int[grid2D->n_segment2d_total * 2];
    ii = 0;
    for(int iLine = 0; iLine < grid2D->lines.size(); iLine++)
    {
        for(int isegment2d = 0;  isegment2d < grid2D->lines[iLine].size(); isegment2d++)
        {
            grid_point0[ii * 2]     = grid2D->lines[iLine][isegment2d].grid_point0[0];
            grid_point0[ii * 2 + 1] = grid2D->lines[iLine][isegment2d].grid_point0[1];
            ii++;
        }
    }

    int *grid_point1 = new int[grid2D->n_segment2d_total * 2];
    ii = 0;
    for(int iLine = 0; iLine < grid2D->lines.size(); iLine++)
    {
        for(int isegment2d = 0;  isegment2d < grid2D->lines[iLine].size(); isegment2d++)
        {
            grid_point1[ii * 2]     = grid2D->lines[iLine][isegment2d].grid_point1[0];
            grid_point1[ii * 2 + 1] = grid2D->lines[iLine][isegment2d].grid_point1[1];
            ii++;
        }
    }

    double *normal = new double[grid2D->n_segment2d_total * 3];
    ii = 0;
    for(int iLine = 0; iLine < grid2D->lines.size(); iLine++)
    {
        for(int isegment2d = 0;  isegment2d < grid2D->lines[iLine].size(); isegment2d++)
        {
            normal[ii * 3]     = grid2D->lines[iLine][isegment2d].normal[0];
            normal[ii * 3 + 1] = grid2D->lines[iLine][isegment2d].normal[1];
            normal[ii * 3 + 2] = grid2D->lines[iLine][isegment2d].normal[2];
            ii++;
        }
    }

    //write n_segment2d
    ndim_segment2d = 1;
    dim1d_segment2d[0] = grid2D->n_line;
    segment2d_dataspace_id = H5Screate_simple(ndim_segment2d, dim1d_segment2d, NULL);
    segment2d_dataset_id = H5Dcreate2(grid_file_id, "n_segment2d", H5T_NATIVE_INT, segment2d_dataspace_id,
                                   H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    H5Dwrite(segment2d_dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, n_segment2d);
    H5Sclose(segment2d_dataspace_id);
    H5Dclose(segment2d_dataset_id);

    //write start_point
    ndim_segment2d = 2;
    dim2d_segment2d[0] = grid2D->n_segment2d_total;
    dim2d_segment2d[1] = 2;
    segment2d_dataspace_id = H5Screate_simple(ndim_segment2d, dim2d_segment2d, NULL);
    segment2d_dataset_id = H5Dcreate2(grid_file_id, "start_point", H5T_NATIVE_DOUBLE, segment2d_dataspace_id,
                                   H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    H5Dwrite(segment2d_dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, start_point);
    H5Sclose(segment2d_dataspace_id);
    H5Dclose(segment2d_dataset_id);
    
    //write end_point
    ndim_segment2d = 2;
    dim2d_segment2d[0] = grid2D->n_segment2d_total;
    dim2d_segment2d[1] = 2;
    segment2d_dataspace_id = H5Screate_simple(ndim_segment2d, dim2d_segment2d, NULL);
    segment2d_dataset_id = H5Dcreate2(grid_file_id, "end_point", H5T_NATIVE_DOUBLE, segment2d_dataspace_id,
                                   H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    H5Dwrite(segment2d_dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, end_point);
    H5Sclose(segment2d_dataspace_id);
    H5Dclose(segment2d_dataset_id);

    //write length
    ndim_segment2d = 1;
    dim1d_segment2d[0] = grid2D->n_segment2d_total;
    segment2d_dataspace_id = H5Screate_simple(ndim_segment2d, dim1d_segment2d, NULL);
    segment2d_dataset_id = H5Dcreate2(grid_file_id, "length", H5T_NATIVE_DOUBLE, segment2d_dataspace_id,
                                   H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    H5Dwrite(segment2d_dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, length);
    H5Sclose(segment2d_dataspace_id);
    H5Dclose(segment2d_dataset_id);

    //write grid_point0
    ndim_segment2d = 2;
    dim2d_segment2d[0] = grid2D->n_segment2d_total;
    dim2d_segment2d[1] = 2;
    segment2d_dataspace_id = H5Screate_simple(ndim_segment2d, dim2d_segment2d, NULL);
    segment2d_dataset_id = H5Dcreate2(grid_file_id, "grid_point0", H5T_NATIVE_INT, segment2d_dataspace_id,
                                   H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    H5Dwrite(segment2d_dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, grid_point0);
    H5Sclose(segment2d_dataspace_id);
    H5Dclose(segment2d_dataset_id);

    //write grid_point1
    ndim_segment2d = 2;
    dim2d_segment2d[0] = grid2D->n_segment2d_total;
    dim2d_segment2d[1] = 2;
    segment2d_dataspace_id = H5Screate_simple(ndim_segment2d, dim2d_segment2d, NULL);
    segment2d_dataset_id = H5Dcreate2(grid_file_id, "grid_point1", H5T_NATIVE_INT, segment2d_dataspace_id,
                                   H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    H5Dwrite(segment2d_dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, grid_point1);
    H5Sclose(segment2d_dataspace_id);
    H5Dclose(segment2d_dataset_id);

    //write normal
    ndim_segment2d = 2;
    dim2d_segment2d[0] = grid2D->n_segment2d_total;
    dim2d_segment2d[1] = 3;
    segment2d_dataspace_id = H5Screate_simple(ndim_segment2d, dim2d_segment2d, NULL);
    segment2d_dataset_id = H5Dcreate2(grid_file_id, "normal", H5T_NATIVE_DOUBLE, segment2d_dataspace_id,
                                   H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    H5Dwrite(segment2d_dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, normal);
    H5Sclose(segment2d_dataspace_id);
    H5Dclose(segment2d_dataset_id);

    grid_status = H5Fclose(grid_file_id);

    delete[] n_segment2d;
    delete[] start_point;
    delete[] end_point;
    delete[] length;
    delete[] grid_point0;
    delete[] grid_point1;
    delete[] normal;

}


void IO2D::readGrid(Grid* grid)
{
    hid_t       grid_dataspace_id;
    hid_t       grid_dataset_id;
    herr_t      grid_status;
    string      grid_dataset_name;

    hsize_t     count[3];              /* size of subset in the file */
    hsize_t     offset[3];             /* subset offset in the file */
    hsize_t     stride[3];
    hsize_t     block[3];

    int ii;
    int grid_ndim = 3;
    hsize_t grid_dims_global[3];

    Grid2D* grid2D = static_cast<Grid2D*>(grid);
    grid_file_name  = "data/grid.h5";
    grid_file_id    = H5Fopen( grid_file_name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);

    grid_dims_global[0] = 1;
    grid_dims_global[1] = grid2D->globalDims_[0];
    grid_dims_global[2] = grid2D->globalDims_[1];


    //=============read grid============================================
    grid_dataset_name = "is_wall";
    grid_dataset_id = H5Dopen2(grid_file_id, grid_dataset_name.c_str(), H5P_DEFAULT);
    grid_status = H5Dread(grid_dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, grid2D->iswall_global_);
    grid_status = H5Dclose(grid_dataset_id);

    grid_dataset_name = "bndr_type";
    grid_dataset_id = H5Dopen2(grid_file_id, grid_dataset_name.c_str(), H5P_DEFAULT);    
    grid_status = H5Dread(grid_dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, grid2D->bndr_global_);
    grid_status = H5Dclose(grid_dataset_id);

    grid_dataset_name = "bndr_val";
    grid_dataset_id = H5Dopen2(grid_file_id, grid_dataset_name.c_str(), H5P_DEFAULT);        
    grid_status = H5Dread(grid_dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, grid2D->bndrVal_global_);
    grid_status = H5Dclose(grid_dataset_id);

    /*
    //===================== read boundary lines ======================================
    hid_t segment2d_type, memtype, point_type, grid_point_type, normal_type;
    hid_t segment2d_dataspace_id, segment2d_dataset_id;
    int ndim_segment2d;
    hsize_t dim1d_segment2d[1];
    hsize_t dim2d_segment2d[2];

    //read n_segment2d
    segment2d_dataset_id = H5Dopen2(grid_file_id, "n_segment2d", H5P_DEFAULT);
    segment2d_dataspace_id = H5Dget_space(segment2d_dataset_id);
    H5Sget_simple_extent_dims(segment2d_dataspace_id, dim1d_segment2d, NULL);
    grid2D->n_line = dim1d_segment2d[0];

    int *n_segment2d = new int[grid2D->n_line];
    H5Dread(segment2d_dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, n_segment2d);
    grid2D->n_segment2d.resize(grid2D->n_line);
    grid2D->n_segment2d_total = 0;
    for(int i = 0; i < grid2D->n_line; i++)
    {
        grid2D->n_segment2d[i] = n_segment2d[i];
        grid2D->n_segment2d_total += n_segment2d[i];
    }
    H5Sclose(segment2d_dataspace_id);
    H5Dclose(segment2d_dataset_id);

    grid2D->lines.resize(grid2D->n_line);
    for(int iLine = 0; iLine < grid2D->n_line; iLine++)
    {
        grid2D->lines[iLine].resize(grid2D->n_segment2d[iLine]);
    }

    //read start_point
    double *start_point = new double[grid2D->n_segment2d_total * 2];
    segment2d_dataset_id = H5Dopen2(grid_file_id, "start_point", H5P_DEFAULT);
    H5Dread(segment2d_dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, start_point);
    H5Dclose(segment2d_dataset_id);

    ii = 0;
    for(int iLine = 0; iLine < grid2D->lines.size(); iLine++)
    {
        for(int isegment2d = 0;  isegment2d < grid2D->lines[iLine].size(); isegment2d++)
        {
            grid2D->lines[iLine][isegment2d].start_point[0] = start_point[ii * 2];
            grid2D->lines[iLine][isegment2d].start_point[1] = start_point[ii * 2 + 1];
            ii++;
        }
    }

    //read end_point
    double *end_point = new double[grid2D->n_segment2d_total * 2];
    segment2d_dataset_id = H5Dopen2(grid_file_id, "end_point", H5P_DEFAULT);
    H5Dread(segment2d_dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, end_point);
    H5Dclose(segment2d_dataset_id);

    ii = 0;
    for(int iLine = 0; iLine < grid2D->lines.size(); iLine++)
    {
        for(int isegment2d = 0;  isegment2d < grid2D->lines[iLine].size(); isegment2d++)
        {
            grid2D->lines[iLine][isegment2d].end_point[0] = end_point[ii * 2];
            grid2D->lines[iLine][isegment2d].end_point[1] = end_point[ii * 2 + 1];
            ii++;
        }
    }

    //read length
    double *length = new double[grid2D->n_segment2d_total];
    segment2d_dataset_id = H5Dopen2(grid_file_id, "length", H5P_DEFAULT);
    H5Dread(segment2d_dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, length);
    H5Dclose(segment2d_dataset_id);

    ii = 0;
    for(int iLine = 0; iLine < grid2D->lines.size(); iLine++)
    {
        for(int isegment2d = 0;  isegment2d < grid2D->lines[iLine].size(); isegment2d++)
        {
            grid2D->lines[iLine][isegment2d].length = length[ii];
            ii++;
        }
    }

    //read grid_point0
    int *grid_point0 = new int[grid2D->n_segment2d_total * 2];
    segment2d_dataset_id = H5Dopen2(grid_file_id, "grid_point0", H5P_DEFAULT);
    H5Dread(segment2d_dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, grid_point0);
    H5Dclose(segment2d_dataset_id);

    ii = 0;
    for(int iLine = 0; iLine < grid2D->lines.size(); iLine++)
    {
        for(int isegment2d = 0;  isegment2d < grid2D->lines[iLine].size(); isegment2d++)
        {
            grid2D->lines[iLine][isegment2d].grid_point0[0] = grid_point0[ii * 2];
            grid2D->lines[iLine][isegment2d].grid_point0[1] = grid_point0[ii * 2 + 1];
            ii++;
        }
    }

    //read grid_point1
    int *grid_point1 = new int[grid2D->n_segment2d_total * 2];
    segment2d_dataset_id = H5Dopen2(grid_file_id, "grid_point1", H5P_DEFAULT);
    H5Dread(segment2d_dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, grid_point1);
    H5Dclose(segment2d_dataset_id);

    ii = 0;
    for(int iLine = 0; iLine < grid2D->lines.size(); iLine++)
    {
        for(int isegment2d = 0;  isegment2d < grid2D->lines[iLine].size(); isegment2d++)
        {
            grid2D->lines[iLine][isegment2d].grid_point1[0] = grid_point1[ii * 2];
            grid2D->lines[iLine][isegment2d].grid_point1[1] = grid_point1[ii * 2 + 1];
            ii++;
        }
    }

    //read normal
    double *normal = new double[grid2D->n_segment2d_total * 3];
    segment2d_dataset_id = H5Dopen2(grid_file_id, "normal", H5P_DEFAULT);
    H5Dread(segment2d_dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, normal);
    H5Dclose(segment2d_dataset_id);

    ii = 0;
    for(int iLine = 0; iLine < grid2D->lines.size(); iLine++)
    {
        for(int isegment2d = 0;  isegment2d < grid2D->lines[iLine].size(); isegment2d++)
        {
            grid2D->lines[iLine][isegment2d].normal[0] = normal[ii * 3];
            grid2D->lines[iLine][isegment2d].normal[1] = normal[ii * 3 + 1];
            grid2D->lines[iLine][isegment2d].normal[2] = normal[ii * 3 + 2];
            ii++;
        }
    }
    grid_status = H5Fclose(grid_file_id);



    delete[] n_segment2d;
    delete[] start_point;
    delete[] end_point;
    delete[] length;
    delete[] grid_point0;
    delete[] grid_point1;
    delete[] normal;

    */
}
